//////////////////////////////////////////////////////////////////////
//                                                                  //
//                                                                  //
//              Copyright (C), 2000  SEIKO EPSON Corp.              //
//              ALL RIGHTS RESERVED                                 //
//                                                                  //
//                                                                  //
//   file name : exp.c                                              //
//                                                                  //
//   Revision history                                               //
//       00/02/10    first release                 M.Igarashi       //
//                                                                  //
//////////////////////////////////////////////////////////////////////


#include <math.h>
#include <smcvals.h>

#define	K0	0.23093347753750233624e-1
#define	K1	0.20202065651286927227886e2
#define	K2	0.1513906799054338915894328e4
#define	S0	0.233184211427481623790295e3
#define	S1	0.4368211662727558498496814e4


// ALGORITHM
// 1.devide x to the multiple of ln2 and the remains 
//      x = xexp * ln2 + g
//       xexp is integer, g is the remains       |g| <= 1/2*ln2
// 2.calculate exponential function of g by next continued fraction
//                         2g       g^2     g^2     g^2     g^2            g^2
//      exp( g ) = 1 +  --------- ------- ------- ------- -------  . . . -------  . . . .
//                       2 - g +   6 +     10 +    14 +    18 +           N  +
// 3.get exponential function of x  
//      exp( x ) = 2^xexp * exp( g )
//

//  the architecture of double floating point
//
//   0 1          11                  31 32                               63 bit
//   -----------------------------------------------------------------------
//  | |   exponent  |                    fraction                           |
//   -----------------------------------------------------------------------
//
//  |               |                   |                                   |
//  |    12 bits          20 bits       |             32 bits               |
//  |            lower word             |            higher word            |
//
//         bit    0         sign bit         (  1 bit  )
//              1 - 11      exponent part    ( 11 bits )
//             12 - 63      fraction part    ( 52 bits )
//
//
 
double exp(double dfX){

	long lX;
	unsigned long ulLx,ulHx;
	double	dfG,dfGt;
	double	dfG2,dfGt2;
	int	iLn;
	double	dfRet;

	
	GETW_L(lX,dfX);		// get low
	GETW_H(ulHx,dfX);		// get high
	
	ulLx = lX&0x7fffffff;					// mask sign

	if( ulLx == 0x00000000 && ulHx==0x00000000 ){
		dfRet = 1.0;						// exp(0.0)=1.0
		return dfRet;
	}
		
	if( ulLx > 0x40863000 ||(ulLx == 0x40863000 && ulHx > 0x0 )){
		if((lX&0x80000000) !=0 ){		// x is too small
			dfRet=0.0;
			return dfRet;
		}else{
			errno = ERANGE;				// value check
			dfRet = P_INF._D;
			return	dfRet;
		}
	}
		
	// calc exp()
	
	dfX*=INVLS;
	iLn=(int)floor(dfX);
	dfG=(dfX-(double)iLn)-0.5;
	dfG2=dfG*dfG;
	dfGt=K0*dfG2+K1;
	dfGt=dfGt*dfG2+K2;
	dfGt*=dfG;
	dfGt2=dfG2+S0;
	dfGt2=dfGt2*dfG2+S1;
	dfRet=dfGt2+dfGt;
	dfRet/=(dfGt2-dfGt);
	dfRet*=RTS;
	dfRet=ldexp(dfRet,iLn);
	
	return dfRet;
}
